{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6ujKLyUcMXW_"
      },
      "outputs": [],
      "source": [
        "%cd /content/drive/MyDrive/大创/CNN-Based-Hyperspectral-Image-Classification-master/"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c7Dv5CvTgAxE"
      },
      "outputs": [],
      "source": [
        "!pip install spectral"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eq1Enc2jbA9u"
      },
      "outputs": [],
      "source": [
        "import keras\n",
        "from keras.layers import Conv2D, Conv3D, Flatten, Dense, Reshape, BatchNormalization\n",
        "from keras.layers import Dropout, Input\n",
        "from keras.models import Model\n",
        "from keras.optimizers import adam_v2\n",
        "from keras.callbacks import ModelCheckpoint\n",
        "from keras.utils import np_utils\n",
        "\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.metrics import confusion_matrix, accuracy_score, classification_report, cohen_kappa_score\n",
        "\n",
        "from operator import truediv\n",
        "\n",
        "from plotly.offline import init_notebook_mode\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import scipy.io as sio\n",
        "import os\n",
        "import spectral\n",
        "\n",
        "os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"0,1\" # 使用编号为0、1号的GPU\n",
        "\n",
        "init_notebook_mode(connected=True)\n",
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_hdsJbzZbHm9"
      },
      "source": [
        "Data Loading"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FwvSdmHDbCkt"
      },
      "outputs": [],
      "source": [
        "## GLOBAL VARIABLES 修改项：数据集的选择、windowsize；添加新数据需在loadData、output_units、reports处相应修改\n",
        "dataset = 'HC'\n",
        "test_ratio = 0.7\n",
        "windowSize = 7\n",
        "def loadData(name):\n",
        "    data_path = os.path.join(os.getcwd(),'data')\n",
        "    if name == 'IP':\n",
        "        data = sio.loadmat(os.path.join(data_path, 'Indian_pines.mat'))['indian_pines']\n",
        "        labels = sio.loadmat(os.path.join(data_path, 'Indian_pines_gt.mat'))['indian_pines_gt']\n",
        "    elif name == 'SA':\n",
        "        data = sio.loadmat(os.path.join(data_path, 'Salinas.mat'))['salinas']\n",
        "        labels = sio.loadmat(os.path.join(data_path, 'Salinas_gt.mat'))['salinas_gt']\n",
        "    elif name == 'PU':\n",
        "        data = sio.loadmat(os.path.join(data_path, 'PaviaU.mat'))['paviaU']\n",
        "        labels = sio.loadmat(os.path.join(data_path, 'PaviaU_gt.mat'))['paviaU_gt']\n",
        "    elif name == 'LK':\n",
        "        data = sio.loadmat(os.path.join(data_path, 'WHU_Hi_LongKou.mat'))['WHU_Hi_LongKou']\n",
        "        labels = sio.loadmat(os.path.join(data_path, 'WHU_Hi_LongKou_gt.mat'))['WHU_Hi_LongKou_gt']\n",
        "    elif name == 'HC':\n",
        "        data = sio.loadmat(os.path.join(data_path, 'WHU_Hi_HanChuan.mat'))['WHU_Hi_HanChuan']\n",
        "        labels = sio.loadmat(os.path.join(data_path, 'WHU_Hi_HanChuan_gt.mat'))['WHU_Hi_HanChuan_gt']\n",
        "    return data, labels\n",
        "def splitTrainTestSet(X, y, testRatio, randomState=345):\n",
        "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=testRatio, random_state=randomState,\n",
        "                                                        stratify=y)\n",
        "    return X_train, X_test, y_train, y_test\n",
        "def applyPCA(X, numComponents=75):\n",
        "    newX = np.reshape(X, (-1, X.shape[2]))\n",
        "    pca = PCA(n_components=numComponents, whiten=True)\n",
        "    newX = pca.fit_transform(newX)\n",
        "    newX = np.reshape(newX, (X.shape[0],X.shape[1], numComponents))\n",
        "    return newX, pca\n",
        "def padWithZeros(X, margin=2):\n",
        "    newX = np.zeros((X.shape[0] + 2 * margin, X.shape[1] + 2* margin, X.shape[2]))\n",
        "    x_offset = margin\n",
        "    y_offset = margin\n",
        "    newX[x_offset:X.shape[0] + x_offset, y_offset:X.shape[1] + y_offset, :] = X\n",
        "    return newX\n",
        "def createImageCubes(X, y, windowSize=5, removeZeroLabels = True):\n",
        "    margin = int((windowSize - 1) / 2)\n",
        "    zeroPaddedX = padWithZeros(X, margin=margin)\n",
        "    # split patches\n",
        "    patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]))\n",
        "    patchesLabels = np.zeros((X.shape[0] * X.shape[1]))\n",
        "    patchIndex = 0\n",
        "    for r in range(margin, zeroPaddedX.shape[0] - margin):\n",
        "        for c in range(margin, zeroPaddedX.shape[1] - margin):\n",
        "            patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   \n",
        "            patchesData[patchIndex, :, :, :] = patch\n",
        "            patchesLabels[patchIndex] = y[r-margin, c-margin]\n",
        "            patchIndex = patchIndex + 1\n",
        "    if removeZeroLabels:\n",
        "        patchesData = patchesData[patchesLabels>0,:,:,:]\n",
        "        patchesLabels = patchesLabels[patchesLabels>0]\n",
        "        patchesLabels -= 1\n",
        "    return patchesData, patchesLabels\n",
        "X, y = loadData(dataset)#加载数据集\n",
        "print('LoadData::data shape='+str(X.shape)+'labels shape='+str(y.shape))#测试代码----------------------\n",
        "\n",
        "\n",
        "K = X.shape[2]\n",
        "K = 30 if dataset == 'IP' else 15\n",
        "X,pca = applyPCA(X,numComponents=K)#PCA降维\n",
        "\n",
        "print(\"PCA::X=\"+str(X.shape)+\"pca=\"+str(pca))#测试代码----------------------\n",
        "\n",
        "X, y = createImageCubes(X, y, windowSize=windowSize)#分小块\n",
        "\n",
        "print(\"createImageCubes::X=\"+str(X.shape)+\"y=\"+str(y.shape))\n",
        "Xtrain, Xtest, ytrain, ytest = splitTrainTestSet(X, y, test_ratio)\n",
        "\n",
        "Xtrain.shape, Xtest.shape, ytrain.shape, ytest.shape\n",
        "print(\"Xtrain:\"+str(Xtrain.shape))\n",
        "print(\"ytrain:\"+str(ytrain.shape))\n",
        "print(\"Xtest:\"+str(Xtest.shape))\n",
        "print(\"ytest:\"+str(ytest.shape))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y-liEO7VboiX"
      },
      "source": [
        "Model and Training"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GQV4S6YDbpKO"
      },
      "outputs": [],
      "source": [
        "Xtrain = Xtrain.reshape(-1, windowSize, windowSize, K, 1)#做3d卷积需要多出一个维度\n",
        "Xtrain.shape\n",
        "\n",
        "ytrain = np_utils.to_categorical(ytrain)\n",
        "ytrain.shape\n",
        "\n",
        "S = windowSize\n",
        "L = K\n",
        "output_units = 9 if (dataset == 'PU' or dataset == 'PC' or dataset =='LK') else 16\n",
        "## input layer\n",
        "input_layer = Input((S, S, L, 1))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cJPmiZcpuzm_"
      },
      "outputs": [],
      "source": [
        "## convolutional layers\n",
        "conv_layer1 = Conv3D(filters=8, kernel_size=(3, 3, 7), padding='same', activation='relu')(input_layer)\n",
        "conv_layer2 = Conv3D(filters=16, kernel_size=(3, 3, 5), activation='relu')(conv_layer1)\n",
        "conv_layer3 = Conv3D(filters=32, kernel_size=(3, 3, 3), activation='relu')(conv_layer2)\n",
        "print(conv_layer3.shape)\n",
        "conv3d_shape = conv_layer3.shape\n",
        "conv_layer3 = Reshape((conv3d_shape[1], conv3d_shape[2], conv3d_shape[3]*conv3d_shape[4]))(conv_layer3)\n",
        "conv_layer4 = Conv2D(filters=64, kernel_size=(3,3), activation='relu')(conv_layer3)\n",
        "\n",
        "flatten_layer = Flatten()(conv_layer4)\n",
        "\n",
        "## fully connected layers\n",
        "dense_layer1 = Dense(units=256, activation='relu')(flatten_layer)\n",
        "dense_layer1 = Dropout(0.4)(dense_layer1)\n",
        "dense_layer2 = Dense(units=128, activation='relu')(dense_layer1)\n",
        "dense_layer2 = Dropout(0.4)(dense_layer2)\n",
        "output_layer = Dense(units=output_units, activation='softmax')(dense_layer2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "aFmfn4-ru8pv"
      },
      "outputs": [],
      "source": [
        "# define the model with input layer and output layer\n",
        "model = Model(inputs=input_layer, outputs=output_layer)\n",
        "model.summary()\n",
        "\n",
        "# compiling the model\n",
        "adam = adam_v2.Adam(learning_rate=0.001, decay=1e-06)\n",
        "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
        "\n",
        "# checkpoint\n",
        "filepath = \"best-model.hdf5\"\n",
        "checkpoint = ModelCheckpoint(filepath, monitor='accuracy', verbose=1, save_best_only=True, mode='max')\n",
        "callbacks_list = [checkpoint]\n",
        "\n",
        "history = model.fit(x=Xtrain, y=ytrain, batch_size=256, epochs=70, callbacks=callbacks_list)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uJrHkEhccEW-"
      },
      "source": [
        "**Validation**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ybxuPukScDOX"
      },
      "outputs": [],
      "source": [
        "# load best weights\n",
        "model.load_weights(\"best-model.hdf5\")\n",
        "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
        "Xtest = Xtest.reshape(-1, windowSize, windowSize, K, 1)\n",
        "Xtest.shape\n",
        "\n",
        "ytest = np_utils.to_categorical(ytest)\n",
        "ytest.shape\n",
        "\n",
        "Y_pred_test = model.predict(Xtest)\n",
        "y_pred_test = np.argmax(Y_pred_test, axis=1)\n",
        "\n",
        "classification = classification_report(np.argmax(ytest, axis=1), y_pred_test, digits=4)\n",
        "print(classification)\n",
        "\n",
        "def AA_andEachClassAccuracy(confusion_matrix):\n",
        "    counter = confusion_matrix.shape[0]\n",
        "    list_diag = np.diag(confusion_matrix)\n",
        "    list_raw_sum = np.sum(confusion_matrix, axis=1)\n",
        "    each_acc = np.nan_to_num(truediv(list_diag, list_raw_sum))\n",
        "    average_acc = np.mean(each_acc)\n",
        "    return each_acc, average_acc\n",
        "\n",
        "def reports (X_test,y_test,name):\n",
        "    #start = time.time()\n",
        "    Y_pred = model.predict(X_test)\n",
        "    y_pred = np.argmax(Y_pred, axis=1)\n",
        "    #end = time.time()\n",
        "    #print(end - start)\n",
        "    if name == 'IP':\n",
        "        target_names = ['Alfalfa', 'Corn-notill', 'Corn-mintill', 'Corn'\n",
        "                        ,'Grass-pasture', 'Grass-trees', 'Grass-pasture-mowed', \n",
        "                        'Hay-windrowed', 'Oats', 'Soybean-notill', 'Soybean-mintill',\n",
        "                        'Soybean-clean', 'Wheat', 'Woods', 'Buildings-Grass-Trees-Drives',\n",
        "                        'Stone-Steel-Towers']\n",
        "    elif name == 'SA':\n",
        "        target_names = ['Brocoli_green_weeds_1','Brocoli_green_weeds_2','Fallow','Fallow_rough_plow','Fallow_smooth',\n",
        "                        'Stubble','Celery','Grapes_untrained','Soil_vinyard_develop','Corn_senesced_green_weeds',\n",
        "                        'Lettuce_romaine_4wk','Lettuce_romaine_5wk','Lettuce_romaine_6wk','Lettuce_romaine_7wk',\n",
        "                        'Vinyard_untrained','Vinyard_vertical_trellis']\n",
        "    elif name == 'PU':\n",
        "        target_names = ['Asphalt','Meadows','Gravel','Trees', 'Painted metal sheets','Bare Soil','Bitumen',\n",
        "                        'Self-Blocking Bricks','Shadows']\n",
        "    elif name == 'LK':\n",
        "        target_names = ['0','1','2','3','4','5','6','7','8']\n",
        "    elif name == 'HC':\n",
        "        target_names = ['0','1','2','3','4','5','6','7','8','9','10','11','12','13','14','15']\n",
        "    \n",
        "    classification = classification_report(np.argmax(y_test, axis=1), y_pred, target_names=target_names)\n",
        "    oa = accuracy_score(np.argmax(y_test, axis=1), y_pred)\n",
        "    confusion = confusion_matrix(np.argmax(y_test, axis=1), y_pred)\n",
        "    each_acc, aa = AA_andEachClassAccuracy(confusion)\n",
        "    kappa = cohen_kappa_score(np.argmax(y_test, axis=1), y_pred)\n",
        "    score = model.evaluate(X_test, y_test, batch_size=32)\n",
        "    Test_Loss =  score[0]*100\n",
        "    Test_accuracy = score[1]*100\n",
        "    \n",
        "    return classification, confusion, Test_Loss, Test_accuracy, oa*100, each_acc*100, aa*100, kappa*100\n",
        "\n",
        "classification, confusion, Test_loss, Test_accuracy, oa, each_acc, aa, kappa = reports(Xtest,ytest,dataset)\n",
        "classification = str(classification)\n",
        "confusion = str(confusion)\n",
        "file_name = \"classification_report.txt\"\n",
        "\n",
        "with open(file_name, 'w') as x_file:\n",
        "    x_file.write('{} Test loss (%)'.format(Test_loss))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Test accuracy (%)'.format(Test_accuracy))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Kappa accuracy (%)'.format(kappa))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Overall accuracy (%)'.format(oa))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Average accuracy (%)'.format(aa))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{}'.format(classification))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{}'.format(confusion))\n",
        "\n",
        "def Patch(data,height_index,width_index):\n",
        "    height_slice = slice(height_index, height_index+PATCH_SIZE)\n",
        "    width_slice = slice(width_index, width_index+PATCH_SIZE)\n",
        "    patch = data[height_slice, width_slice, :]\n",
        "    \n",
        "    return patch\n",
        "\n",
        "# load the original image\n",
        "X, y = loadData(dataset)\n",
        "\n",
        "height = y.shape[0]\n",
        "width = y.shape[1]\n",
        "PATCH_SIZE = windowSize\n",
        "numComponents = K\n",
        "\n",
        "X,pca = applyPCA(X, numComponents=numComponents)\n",
        "\n",
        "X = padWithZeros(X, PATCH_SIZE//2)\n",
        "\n",
        "# calculate the predicted image\n",
        "outputs = np.zeros((height,width))\n",
        "for i in range(height):\n",
        "    for j in range(width):\n",
        "        target = int(y[i,j])\n",
        "        if target == 0 :\n",
        "            continue\n",
        "        else :\n",
        "            image_patch=Patch(X,i,j)\n",
        "            X_test_image = image_patch.reshape(1,image_patch.shape[0],image_patch.shape[1], image_patch.shape[2], 1).astype('float32')                                   \n",
        "            prediction = (model.predict(X_test_image))\n",
        "            prediction = np.argmax(prediction, axis=1)\n",
        "            outputs[i][j] = prediction+1\n",
        "\n",
        "ground_truth = spectral.imshow(classes = y,figsize =(7,7))\n",
        "\n",
        "predict_image = spectral.imshow(classes = outputs.astype(int),figsize =(7,7))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XosIzeHbcSEO"
      },
      "source": [
        "**save**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OuIX2k6_cRtt"
      },
      "outputs": [],
      "source": [
        "spectral.save_rgb(\"predictions.jpg\", outputs.astype(int), colors=spectral.spy_colors)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "Untitled1.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "nbformat": 4,
  "nbformat_minor": 0
}